import random
from matplotlib import colors, interactive, scale
import numpy
import math
import matplotlib.pyplot as plt
import copy
#interval
a = 0
b = 10
#parameter
c = 1.0
d = 2.0
e = 3.0
f = 4.0
#选取的采样点数量，最总的插值函数的最高次项 = N -1 需要计算 N 个方程组
N = 6
# 0 -> 等距插值 1 -> 随机插值
israndom = 1

interval = (a, b)
#绘图的粒度
drawscaleX = numpy.arange((interval[0]), (interval[1]), 0.01)
drawPoint  = []
#Vandermonde
VandermondeCoeffientList = []
#牛顿插值法
DifferenceQuotientList= []
linerinter = lambda x0,x1,x:((((x-x1)/(x0-x1))*getFuctionValue(x0)) + (((x-x0)/(x1-x0))*getFuctionValue(x1)))
#Cubic Hermite
deltax = 0.00000001
alpha  = lambda x0,x1,x:(math.pow(((x-x1)/(x0-x1)),2)*(1+2*((x-x0)/(x1-x0)))*getFuctionValue(x0))
beta   = lambda x0,x1,x:(math.pow(((x-x1)/(x0-x1)),2))*(x-x0)*derivative(x0)

step = (interval[1] - interval[0]) / N
EffectNumber = 3

pointList = []

def getFuctionValue(x):
    value = (c * math.sin(d * x)) + (e * math.cos(f * x))
    return value

def derivative(x):
    global deltax
    return ((getFuctionValue(x+deltax) - getFuctionValue(x))/deltax)

def initList():
    global drawPoint
    if israndom:
        pointList.append(interval[0])
        pointList.append(interval[1])
        for i in range(0,N - 2):
            x = round((random.uniform(interval[0], interval[1])), EffectNumber)
            pointList.append(x)
        pointList.sort()
    else:
        step = (interval[1] - interval[0])/(N - 1)
        index = interval[0]
        for i in range(0,N - 1):
            pointList.append(round(index,EffectNumber))
            index += step
        pointList.append(interval[1])
    drawPoint = drawscaleX.tolist()
    for point in pointList:
        drawPoint.append(point)
    drawPoint.sort()
    return

def OmigaX(x, index):
    if index == 0:
        return 1
    result = 1
    for i in range(0,index):
        result *= (x - pointList[i])
    return result

def getVandermondeCoefficient(Vlist:list):
    A = numpy.ndarray(shape=(N,N))
    for i in range(0,N):
        for k in range(0,N):
            A[i][k] = math.pow(pointList[i], k)
    B = numpy.ndarray(shape=(N))
    for i in range(0,N):
        B[i] = getFuctionValue(pointList[i])
    result = numpy.linalg.solve(A,B)
    for i in range(0,N):
        Vlist.append(result[i])
    return

def Vandermonde_Interpolation(x):
    result = 0
    for i in range(0,N):
        result += VandermondeCoeffientList[i] * math.pow(x,i)
    return result

def Lagrange_Interpolation(x):
    result = 0
    for k in range(0, N):  # k : 拉格朗日函数的第k项
        numerator = 1  #分子
        denominator = 1  # 分母
        for i in range(0, N):
            if i == k:
                pass
            else:
                numerator   = numerator   * (x            - pointList[i])
                denominator = denominator * (pointList[k] - pointList[i])
        result += (numerator / denominator) * getFuctionValue(pointList[k])
    return result


def getDQList(DQList:list):
    DQTabel = []
    getDQ = lambda column,raw:(DQTabel[column -1 ][raw + 1] - DQTabel[column - 1][raw])/(pointList[column] - pointList[0])
    for i in range(0,N):
        DQTabelColumn = []
        for k in range(0,N-i):
            if(i == 0):
                DQTabelColumn.append(getFuctionValue(pointList[k]))
            else:
                DQTabelColumn.append(getDQ(i,k))
        DQTabel.append(DQTabelColumn)
    for i in range(0,N):
        DQList.append(DQTabel[i][0])
    return

#对于 N 个插值点，可以计算其 N - 1 阶差商，先计算差商表。
def Newton_Interpolation(x):
    result = 0
    for i in range(0,N):
        result += DifferenceQuotientList[i]*OmigaX(x,i)
    return result

def Liner_Interpolation(x):
    index = 0
    if(x <= pointList[index]):
        index = 0
    elif(x >= pointList[N-1]):
        index = N-2
    else:
        while(not (x > pointList[index] and x <= pointList[index + 1])):
            index += 1
    return linerinter(pointList[index],pointList[index+1],x)

def cubicHermite_Interprolation(x):
    index = 0
    if(x <= pointList[index]):
        index = 0
    elif(x >= pointList[N-1]):
        index = N-2
    else:
        while(not (x > pointList[index] and x <= pointList[index + 1])):
            index += 1
    return  alpha(pointList[index],pointList[index+1],x) + \
            alpha(pointList[index+1],pointList[index],x) + \
            beta(pointList[index],pointList[index+1],x)  + \
            beta(pointList[index+1],pointList[index],x) 

def draw(funcname, labelstr , color):
    resultList = []
    for point in drawPoint :
        resultList.append(funcname(point))
    plt.plot(drawPoint, resultList, color , label = labelstr)
    samplingSpots = []
    for point in pointList:
        samplingSpots.append(funcname(point))
    plt.scatter(pointList, samplingSpots , s=25, c = color)

if __name__ == "__main__":
    initList()
    getVandermondeCoefficient(VandermondeCoeffientList)
    getDQList(DifferenceQuotientList)
    draw(Vandermonde_Interpolation,"Vandermonde",'r')
    draw(Newton_Interpolation,"Newton",'b')
    draw(Lagrange_Interpolation,"Lagrange",'y')
    draw(Liner_Interpolation,"Liner",'k')
    draw(cubicHermite_Interprolation,"cubicHermite",'c')
    draw(getFuctionValue,"origin",'g')
    plt.legend()
    plt.show()